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ABSTRACT 

Wavelet-based  distributed  data  processing  holds  much  promise  for 
sensor  networks;  however,  irregular  sensor  node  placement  pre¬ 
cludes  the  direct  application  of  standard  wavelet  techniques.  In 
this  paper,  we  develop  a  new  distributed  wavelet  transform  based 
on  lifting  that  takes  into  account  irregular  sampling  and  provides  a 
piecewise-planar  multiresolution  representation  of  the  sensed  data. 
We  develop  the  transform  theory;  outline  how  to  implement  it  in 
a  multi-hop,  wireless  sensor  network;  and  illustrate  with  several 
simulations.  The  new  transform  performs  on  par  with  conventional 
wavelet  methods  in  a  head-to-head  comparison  on  a  regular  grid  of 
sensor  nodes. 

1.  INTRODUCTION 

Wireless  sensor  networks  have  emerged  as  an  important  applica¬ 
tion  area  for  distributed  signal  processing.  Sensor  networks  consist 
of  nodes  that  sense  phenomena  of  interest,  process  the  measure¬ 
ments,  and  share  data  via  a  wireless,  multi-hop  routing  network. 
Nodes  have  limited  on-board  power  supplies,  and  since  commu¬ 
nication  power  consumption  typically  dominates  over  processing 
power  by  orders  of  magnitude,  intelligent,  in-network  signal  pro¬ 
cessing  is  necessary  to  reduce  the  amount  of  transmitted  data. 
Whenever  possible,  transmissions  outside  the  network  should  take 
the  form  of  summarized  results  and  conclusions  rather  than  raw 
data.  Such  processing  must  be  both  distributed  —  not  requiring 
all  data  in  a  central  location  —  and  localized  —  requiring  access 
only  to  data  in  a  node’s  immediate  vicinity. 

The  restrictions  on  signal  processing  algorithms  for  sensor 
networks  are  considerably  complicated  by  the  irregular  node 
placement  typical  of  real-world  deployments.  While  most  tradi¬ 
tional  regular-grid  signal  processing  techniques  do  not  translate 
directly  to  this  setting,  much  of  the  literature  on  distributed  sig¬ 
nal  processing  for  sensor  networks  has  nonetheless  assumed  reg¬ 
ular  sampling  grids,  a  fact  highlighted  in  [1],  To  help  rectify  this 
discrepancy,  we  propose  in  this  paper  what  is  to  our  knowledge 
the  first  distributed,  two  dimensional  (2-D),  irregular-grid  wavelet 
transform  for  sensor  networks  capable  of  multiscale,  piecewise- 
planar  approximation  of  node  measurements  (two  wavelet  vanish¬ 
ing  moments).  We  provide  a  detailed  treatment  not  only  of  the 
transform  theory  (based  on  wavelet  lifting  [2])  but  also  its  imple¬ 
mentation  issues,  developing  along  the  way  a  new  distributed  tri¬ 
angulation  protocol  that  extends  the  work  of  [3]. 

Section  2  discusses  related  work  in  wavelet  processing  for  sen¬ 
sor  networks.  Section  3  overviews  the  lifting  approach  that  we 
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exploit  in  the  later  sections.  Section  4  develops  a  multiscale  struc¬ 
ture  on  the  sensor  nodes;  Section  5  derives  the  requisite  filters; 
and  Section  6  deals  with  the  iteration  of  the  transform.  Section 
7  develops  a  new  technique  for  distributed  network  triangulation. 
Section  8  demonstrates  the  approximation  power  of  our  transform 
coefficients  with  two  distributed  compression  examples.  Section  9 
concludes  and  reviews  our  ongoing  work. 

2.  RELATED  WORK 

Multiscale  algorithms  are  not  difficult  to  motivate  for  sensor  net¬ 
work  applications,  since  the  laws  of  physics  often  induce  in  the 
measured  data  a  natural  multiscale  structure  that  can  guide  the 
scope  and  extent  of  in-network  signal  processing  and  communi¬ 
cation.  In  particular,  as  sensors  become  more  distant  from  each 
other,  the  spatial  correlations  between  their  measurements  will  de¬ 
cay  rapidly.  This  suggests  local  processing  at  fine  scales  between 
neighboring  nodes  and  global  processing  at  coarse  scales  between 
more  far-flung  nodes. 

DIMENSIONS  [4]  uses  an  in-network  wavelet  transform  to  fa¬ 
cilitate  querying  and  storage  of  sensor  network  measurements,  but 
it  assumes  a  regular-grid  placement  of  nodes.  The  same  assump¬ 
tion  is  shared  by  the  wavelet-based  Wisden  system  [5]  for  struc¬ 
tural  monitoring.  Similarly,  [6]  proposes  separable  application  of 

1- D  regular-grid  wavelet  transforms  to  solve  the  2-D  sensor  broad¬ 
cast  problem.  The  lifting-based,  regular-grid  distributed  wavelet 
transform  in  [7]  is  similar  in  spirit  to  the  one  proposed  here.  It 
employs  a  1-D  wavelet  decomposition  along  a  path  through  the 

2- D  measurement  field;  however,  no  method  for  determining  the 
optimal  path  is  given.  While  the  technique  could  be  extended 
to  use  irregular-grid  1-D  wavelets,  such  an  approach  is  not  capa¬ 
ble  of  fully  capturing  the  higher-dimensional  spatial  dependencies 
among  the  measurements.  Similar  conclusions  apply  to  the  1-D 
Haar  protocol  described  in  [8]. 

The  work  in  [9],  which  this  paper  extends,  provides  an 
irregular-grid,  fully  2-D,  distributed  wavelet  transform  for  sensor 
networks,  based  on  piecewise-constant  multiscale  approximation 
and  multiscale  routing  structures.  In  this  paper,  we  develop  a 
distributed  lifting  transform  capable  of  piecewise-planar  approx¬ 
imation  and  requiring  no  a  priori  multiscale  network  structure. 
There  has  been  significant  prior  treatment  of  centralized  irregular- 
grid  lifting  in  both  the  computer  graphics  and  statistical  estimation 
communities  (see  [10, 1 1]  and  the  references  therein);  we  base  our 
distributed  scheme  here  on  a  technique  suggested  in  [1 1]. 

3.  WAVELET  LIFTING  ON  IRREGULAR  GRIDS 

Wavelet  lifting  [2]  replaces  the  measurement  at  each  sensor  net¬ 
work  node  with  a  wavelet  coefficient  representing  the  inner  prod- 
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uct  of  the  measured  field  with  a  wavelet  function  in  a  basis  or  frame 
expansion.  The  process  is  multiscale,  starting  from  the  original 
node  measurements  at  some  finest  scale  J  and  iterating  to  a  final, 
coarsest  scale  0.  Denote  the  complete  set  of  nodes  by  E(  J)  and 
the  complete  measurements  by  sj\  these  comprise  the  finest-scale 
set  of  scaling  coefficients.  We  assume  that  each  node  knows  its 
location  in  space;  providing  this  information  is  an  area  of  active 
research  [12], 

Adopting  with  slight  modification  the  notation  of  [11]  and 
starting  with  j  =  J  —  1,  we  say  that  each  scale  j  of  the  transform 
first  splits  the  set  of  nodes  E(j  +  1)  into  a  set  of  ff=E(j)  “even” 
coefficients  located  at  nodes  E(j)  and  a  set  of  jfO(j)  “odd”  co¬ 
efficients  located  at  nodes  O(j).1  This  divides  the  set  of  scaling 
coefficients  into  two  new  sets  •Sj+i,bo)  and  .Sj+i.oO')-  The  co¬ 
efficients  in  Sj+i,o(j)  give  rise  t°  the  scale- j  wavelet  coefficients 
We  say  that  these  are  predicted  using  Sj+itE(j)  and  re¬ 
fer  to  the  nodes  in  0(j )  as  predicted  nodes. The  coefficients  in 
sj+i,E(j)  are  then  updated  to  SjtE(j)  using  d:ho(j)  in  order  to  pre¬ 
serve  the  weighted  average  of  the  scaling  coefficients;  the  nodes  in 
E(j)  are  known  as  updated  nodes. 

The  scale- j  set  of  scaling  coefficients  SjtEy)  replaces  the  even 
scaling  coefficients  from  scale  j  - f  1,  and  the  transform  then  iter¬ 
ates  over  s jjE(j)  at  coarser  scale  j  —  1,  continuing  until  a  root 
set  of  scaling  coefficients  resides  at  the  nodes  in  E( 0).  At  that 
point,  the  set  of  terminal  scaling  coefficients  so  and  the  entire  set 
of  wavelet  coefficients  {d?}jg{o,...,  J-i}  number  the  same  as  the 
original  field  measurements  at  E(J)  and  completely  describe  them 
—  that  is,  the  entire  process  is  invertible. 

The  key  issues  in  any  lifting  decomposition  are  (1)  determin¬ 
ing  which  values  are  to  be  decimated  at  each  scale  (i.e.,  partition¬ 
ing  nodes  into  E(j)  and  O(j)),  and  (2)  defining  filters  to  calculate 
the  scaling  and  wavelet  coefficients  at  each  scale.  In  a  regular-grid 
setting,  the  grid  structure  guides  both  these  choices  —  decima¬ 
tion  to  a  coarser  grid  is  trivial,  and  the  same  set  of  scaling  and 
wavelet  filters  can  be  applied  at  each  grid  point.  With  grid  irreg¬ 
ularity,  choosing  which  nodes  to  decimate  becomes  more  compli¬ 
cated,  and  the  lack  of  a  predictable  spatial  neighbor  distribution 
necessitates  adapting  different  filter  coefficients  for  each  node. 

The  technique  in  [9]  exploits  the  simplicity  of  piecewise- 
constant  multiscale  approximation  and  leverages  hierarchical  rout¬ 
ing  structures  (see  [13]  for  an  example)  to  sidestep  these  two  is¬ 
sues.  The  routing  groups  nodes  into  clusters,  and  the  transform 
assigns  within  each  cluster  a  scaling  coefficient  (average)  to  an 
elected  clusterhead  node  while  assigning  wavelet  coefficients  (dif¬ 
ferences  from  the  average)  to  the  other  nodes.  Clusterheads  partic¬ 
ipate  in  the  next  transform  level.  The  routing  hierarchy  thus  guides 
decimation,  and  simple  averaging  avoids  the  need  for  filters  that 
depend  on  spatial  neighbor  distributions.  Unfortunately,  when  the 
data  are  smoother  than  piecewise-constant,  this  approximation  will 
not  completely  sparsify  the  wavelet  coefficients,  degrading  perfor¬ 
mance  in  applications  such  as  compression  and  de-noising.  To  en¬ 
able  higher-order  multiscale  approximation  and  to  free  ourselves 
from  dependency  on  multiscale  routing  structures,  we  must  meet 
both  the  challenges  of  decimation  and  local  filter  design  head-on. 

4.  SENSOR  NODE  DECIMATION 

At  each  scale  j  we  separate  predicted  nodes  from  the  previous 
scale  into  predicted  nodes  O(j)  and  updated  nodes  E(j),  where 

1  #  denotes  here  the  number  of  elements  in  a  set. 


E(j  +  1)  =  E(j)  U  O(j).  Naturally,  0(j)  and  E(J)  must  be 
disjoint.  Additionally,  we  wish  the  decimated  set  O(j)  to  be  as 
large  as  possible  at  each  scale  in  order  to  have  as  few  scales  as 
possible.  To  accomplish  this,  we  construct  a  mesh  to  guide  the 
decimation  —  in  particular  (but  without  loss  of  generality)  a  De¬ 
launay  triangulation  (DT)  of  the  nodes  at  scale  j  [10,  11].  DT’s 
can  be  built  in  a  distributed  fashion,  as  detailed  in  Section  7.  The 
triangulation  is  calculated  once  (at  the  finest  level  of  the  transform) 
and  then  locally  updated  as  points  are  removed  at  each  scale.  To 
begin,  each  node  need  only  discover  its  location  and  share  this  in¬ 
formation  among  nodes  within  its  communication  range  (see  [12] 
and  the  references  therein). 

Using  the  mesh,  we  first  define  the  neighbors  of  a  node  at  scale 
j  as  those  nodes  with  which  it  shares  an  edge  in  the  scale- j  mesh. 
We  can  then  designate  which  points  to  remove  at  each  transform 
level  by  appealing  to  a  notion  of  scale  in  the  irregular  grid.  The 
mesh  allows  us  to  assign  an  area  of  support  to  each  node  —  a 
computation  detailed  in  Section  5.3.  While  in  a  regular  grid  each 
node  has  the  same  support  area,  each  node  in  an  irregular  grid  has 
a  different  area  and,  thus,  a  unique  scale.  Since  multiscale  anal¬ 
ysis  provides  increasingly  broad  views  of  a  measurement  field  at 
each  scale,  it  follows  naturally  that  we  should  remove  the  smallest- 
scale  (-area)  points  at  each  transform  level.  Recall,  however,  that 
a  node  predicted  at  level  j  cannot  have  predicted  neighbors  at  that 
level.  To  avoid  violating  these  conditions,  we  propose  the  follow¬ 
ing  greedy  solution  to  decimating  nodes: 

1.  Mark  all  nodes  on  the  mesh  boundary  as  updated,  placing 
them  in  E(j). 

2.  Mark  the  node  with  the  smallest  area  among  unmarked 
nodes  as  predicted,  placing  it  in  O(j). 

3.  Mark  the  predicted  node’s  neighbors  in  the  mesh  as  up¬ 
dated,  placing  them  in  E(j). 

4.  Return  to  Step  2  while  unmarked  nodes  remain. 

Step  1  ensures  the  stability  of  the  predict  stage  (discussed 
further  in  Section  5.2)  and  guarantees  that  boundary  nodes  form 
the  coarsest  set  .E(O).  While  this  protocol  is  essentially  a  cen¬ 
tralized  one,  it  provides  an  elegant  solution  to  decimation  with 
relatively  low  overhead.  One  could  imagine  a  message-passing 
scheme  where  nodes  eventually  agree  on  which  one  has  the  next 
smallest  area;  however,  the  associated  communication  overhead 
could  become  prohibitive  as  the  network  grows.  The  triangulation 
algorithm  is  deterministic  given  node  locations,  and  so  we  can  eas¬ 
ily  compute  the  decimation  order  outside  the  network  after  gather¬ 
ing  the  nodes’  self-localized  positions  and  then  inform  individual 
nodes  of  the  transform  level  at  which  they  are  decimated. 

Figure  1  illustrates  two  successive  levels  of  decimation.  Pre¬ 
dicted  nodes  at  scale  j  are  marked  with  dark  circles  in  Figure  1(a). 
Following  removal  of  these  points  and  local  re-triangulation,  a  new 
set  of  predicted  nodes  is  designated  at  scale  j  —  1  in  Figure  1(b). 
Experimentally,  this  technique  removes  approximately  25%  of  the 
nodes  per  scale. 

5.  DISTRIBUTED  FILTER  DESIGN 

Once  the  sets  O(j)  and  E(j)  have  been  designated,  we  must  com¬ 
pute  the  scale-y  wavelet  coefficients  for  the  nodes  in  O(j)  and 
scaling  coefficients  for  the  nodes  in  E(j).  This  entails  assigning 
predict  and  update  filters  to  operate  on  the  coefficients  at  neigh¬ 
boring  nodes.  When  the  node  locations  remain  fixed,  both  the  dec¬ 
imation  order  and  filter  coefficients  need  only  be  computed  once. 


Fig.  1.  Triangulated  mesh  of  sensor  nodes  (a)  at  scale  j  and  (b)  at 
coarser  scale  ji  —  If*  denotes  vertices  to  be  decimated). 


While  we  suggest  that  the  decimation  structure  be  centrally  com¬ 
puted  and  disseminated  to  the  nodes  as  outlined  in  Section  4,  the 
filter  coefficients  can  be  efficiently  computed  in  a  distributed  fash¬ 
ion  and  stored  locally  at  the  nodes  for  future  reference. 

5.1.  Refinement  relations 

We  first  specify  how  the  wavelet  and  scaling  functions  are  related 
across  scales  in  our  piecewise-planar  biorthogonal  wavelet  trans¬ 
form,  which  is  inspired  by  Section  7.2.3  of  [11].  Label  the  row 
vector  of  scaling  functions  to  be  generated  at  scale  j  for  nodes 
in  E(j)  as  &jtE(j),  and  let  the  scaling  coefficients  SjtE(j)  form 
a  column  vector.  Similarly,  let  'Lj,o(j)  represent  the  row  vector 
of  scale-/  wavelet  functions  at  nodes  in  O(j)  with  an  associated 
column  vector  dj^o(j)  of  wavelet  coefficients.  To  relate  the  scal¬ 
ing  functions  from  the  finer  scale  /  +  1  to  the  scaling  and  wavelet 
functions  at  scale  /,  define  an  #0(j)  x  #E(j)  predict  filter  Pj 
and  an  #E(j)  x  #0(j)  update  filter  Uj  to  operate  as  follows 


regress  a  plane  through  its  neighbors  in  the  scale-/  mesh.  De¬ 
note  these  neighbors  as  A f(np)  £  E(j).  Per  (2),  the  value  of 
the  wavelet  coefficient  is  the  difference  between  the  scale- (j  +  1) 
scaling  coefficient  at  np  and  the  value  predicted  by  the  plane  at  the 
location  of  nv.  The  1  x  f(np)  vector  p,  of  predict  coeffi- 

— jinp 

cients2  at  np  is  given  by 

p  ={l,x(np),y(np)](X'X)-1X',  (3) 

— J  >  ,Lp 

where  x(-)  and  y(-)  give  the  x  and  y  coordinates  of  a  node  and 
the  #Af(np)  x  3  matrix  .Y  =  [1,  x(Af(np)),y(Af(np))\  (where 
1  is  a  column  vector  of  ones).  Properly  sorted,  the  coefficients 
in  p .  fill  sensor  np  s  row  in  Pj  at  the  positions  corresponding 

— 3inp 

to  neighbors  Af(np).  The  remainder  of  the  elements  in  the  row, 
corresponding  to  more  distant  nodes,  are  zero,  and  thus  the  entire 
process  is  local.  Figure  2(a)  depicts  the  one-time  local  data  flow 
at  scale  /  during  computation  of  predict  filter  coefficients,  with 
neighbors  of  predicted  nodes  transmitting  their  ( x ,  y)  coordinates 
for  use  in  (3). 

This  predict  scheme  illuminates  why  care  is  taken  in  Section 
4  to  exclude  boundary  nodes  when  selecting  nodes  for  prediction. 
When  the  angular  distribution  in  polar  coordinates  of  neighbors 
around  a  predicted  node  does  not  cover  most  of  the  range  [0,  27r), 
planar  regression  becomes  extrapolation  in  the  region  with  no  data 
points,  leading  to  numerical  instability.  Several  intricate,  central¬ 
ized  techniques  are  mentioned  in  [11]  to  deal  with  this  instabil¬ 
ity.  To  counter  this  problem  in  a  distributed  fashion,  we  simply 
choose  never  to  decimate  nodes  on  the  boundary  of  the  finest-level 
mesh,  since  these  are  precisely  the  nodes  that  do  not  have  neigh¬ 
bors  throughout  the  surrounding  [0,  2n)  space. 


-  $j+l,E(j)  +  $j+i ,o(j)Pj  (1) 

=  ®j+i,oU)  ~  ,E(j)Uj. 

Similarly,  Pj  and  Uj  relate  the  scaling  coefficients  from  scale  /  + 1 
to  the  scaling/wavelet  coefficients  at  coarser  scale  /  as 

dj,0(j)  =  SJ  +  l,0(j)  —  PjSjj. l,E(j) 

Sj,E(j)  =  sj+l,E(j)  +  Ujdjy0(j)- 

Without  careful  design  the  transform  may  be  far  from  orthogo¬ 
nal,  resulting  in  coefficients  whose  impact  on  reconstruction  in  the 
spatial  domain  is  poorly  related  to  their  magnitudes  in  the  wavelet 
domain.  Entries  in  Pj  and  Uj  must  therefore  be  chosen  so  that  the 
transform  is  as  close  to  orthogonal  as  possible.  The  filters  also  con¬ 
trol  how  distributed  and  local  the  transform  will  be.  An  element  in 
a  row  of  Pj  or  Uj  must  be  nonzero  only  if  the  node  corresponding 
to  that  element’s  column  is  a  neighbor  of  the  node  corresponding 
to  the  element’s  row  in  the  scale-/  mesh.  Moreover,  data  gath¬ 
ered  in  calculating  these  predict  and  update  filter  coefficients  at 
each  node  must  come  only  from  neighbors  of  each  node  at  scale  /. 
With  these  requirements  in  mind,  we  describe  how  to  locally  fill  in 
the  entries  of  P,  and  Uj  and  discuss  the  associated  communication 
traffic. 

5.2.  Predict  filters 

The  predict  filter  has  the  most  straightforward  design.  We  desire 
the  transform  to  have  two  vanishing  moments;  that  is,  if  the  data 
are  planar,  then  the  wavelet  coefficients  should  be  zero.  To  com¬ 
pute  the  wavelet  coefficient  at  a  predicted  point  np  £  0(j),  we 


5.3.  Scaling  function  support  areas 

Computing  the  update  coefficients,  to  follow  shortly,  requires 
maintaining  integrals  of  the  scaling  functions  at  scale  /  in  terms  of 
integrals  of  the  scaling  functions  at  scale  j  +  1.  If  we  consider  the 
scaling  functions  at  the  finest  scale  to  be  indicator  functions  over 
pseudo-Voronoi  regions  surrounding  the  nodes,  then  this  amounts 
to  calculating  and  updating  the  scaling  function  support  areas.  The 
first  set  of  integrals  is  extracted  from  the  finest-scale  mesh,  with 
each  node  assigning  to  itself  1/3  the  area  of  each  triangle  to  which 
it  belongs  [10], 

Denote  by  AjtUu  the  integral  of  scaling  function  (f)j,nu  £ 
&j,E(j)  at  updated  node  nu  £  Ej.  Integrating  the  first  equation  of 
(1)  at  nu  yields 

Aj,nu  —  Aj+l,nu  +  y  '  Pj  (nu)Aj+l,nk,  (4) 

* — JULk 

nkeU(nu) 

where  A f(nu)  describes  the  neighbors  of  nu  whose  wavelet  coef¬ 
ficients  are  predicted  at  scale  /;  p.  ( nu )  gives  the  predict  co¬ 
efficient  node  nk  applies  to  the  scale-(/  +  1)  scaling  coefficient 
of  node  nu.  Essentially,  to  compute  the  integral  of  the  scale-/ 
scaling  function  for  an  updated  node  n„,  we  add  the  integral  of 
its  seal e-(/  +  1)  scaling  function  to  a  weighted  sum  of  the  scale- 
(/  +  1)  scaling  function  integrals  of  its  predicted  neighbors.  The 
prediction  coefficients  for  nu  at  those  neighbors  provide  the  sum 
weights.  Figure  2(b)  illustrates  the  one-time  communication  traffic 

Underlining  is  used  here  to  distinguish  locally  computed  filter  coeffi¬ 
cient  vectors. 


Fig.  3.  Communication  at  scale  j  between  predicted  (•)  and  up¬ 
dated  (o)  nodes  repeated  during  each  wavelet  transform.  Messages 
are  labelled  between  a  predicted  node  np  and  an  updated  neighbor 
nu.  (a)  nu  sends  its  scale-(jf  +  1)  scaling  coefficient  Sj+i,n„  to 
each  predicted  neighbor;  (b)  np  sends  its  scale- j  wavelet  coeffi¬ 
cient  dj:np  value  to  each  updated  neighbor. 


where  d:hJ\f(n\  is  a  column  vector  of  wavelet  coefficients.  Fig¬ 
ure  3(b)  illustrates  the  associated  traffic  pattern,  with  each  updated 
node  passing  its  scale- j  wavelet  coefficient  to  each  of  its  neigh¬ 
bors  to  enable  their  scaling  values  calculations. 

Again,  no  centralized  administration  is  necessary  to  manage 
this  process,  apart  from  all  nodes  in  the  network  agreeing  to  be¬ 
gin  a  new  transform  calculation  —  a  process  that  can  proceed  at 
a  time-scale  appropriate  for  the  measured  phenomenon  and/or  op¬ 
erator  interest.  When  the  measured  field  is  transformed  to  scale 
j,  each  predicted  node  np  merely  waits  for  nodes  in  A f{np)  to 
send  their  scale-(y  +  1)  scaling  coefficients  before  computing  its 
scale-j  wavelet  coefficient.  Similarly,  each  Jin  waits  for  the  scale- 
j  wavelet  coefficients  from  its  predicted  neighbors  before  comput¬ 
ing  its  new  scaling  coefficient.  The  transform  proceeds  through 
scales  in  a  completely  asynchronous  and  distributed  fashion,  never 
progressing  to  the  next,  coarser  scale  until  all  computations  at  the 
previous,  finer  scale  are  complete. 

7.  DISTRIBUTED  TRIANGULATION 

As  detailed  in  Section  3,  a  mesh  is  required  at  each  scale  of  the 
transform  to  provide  connectivity  for  the  predict  and  update  calcu¬ 
lations,  and  areas  extracted  from  the  mesh  at  the  starting  scale  are 
required  to  initialize  the  scaling  function  integral  re-calculation  of 
Section  5.3.  For  the  lifting  transform  to  be  truly  distributed,  con¬ 
struction  and  maintenance  of  the  mesh  must  be  distributed. 

With  some  modifications,  the  algorithm  from  [3]  suffices  to 
compute  a  triangulated  mesh.  The  algorithm  builds  a  Planar  Local¬ 
ized  Delaunay  triangulation  (PLDel).  Given  a  node  communica¬ 
tion  radius,  the  resulting  PLDel  contains  the  edges  in  a  centralized 
DT  whose  lengths  do  not  exceed  the  communication  radius.  First, 
nodes  broadcast  their  locations,  and  each  node  collects  the  loca¬ 
tions  of  neighbors  within  its  communication  radius.  Each  node 
then  locally  computes  a  Delaunay  triangulation  based  on  its  col¬ 
lected  points  and  broadcasts  the  associated  triangles  to  its  neigh¬ 
bors  within  communication  range.  Once  a  node  has  received  these 
triangle  messages,  it  consolidates  message  triangles  with  locally 
computed  triangles  to  obtain  a  final  set  of  valid  Delaunay  triangles 
with  vertices  at  that  node.  The  triangle  edges  with  the  node  as  an 
endpoint,  along  with  any  non-included  Gabriel  edges,3  form  the 

3  A  Gabriel  edge  is  one  for  which  the  circle  whose  diameter  coincides 


Fig.  4.  Test  data  fields:  (a)  noisy  quadratic  bump  with  discontinu¬ 
ity,  and  (b)  Gaussian  bumps  on  a  smooth,  quadratic  field. 


set  of  PLDel  edges  rooted  at  that  node.  Edge  discovery  is  guaran¬ 
teed  to  be  two-way,  meaning  that  each  edge  is  discovered  by  both 
of  its  vertices. 

We  require  more  than  just  connectivity,  however.  A  node  must 
also  know  about  the  triangles  rooted  at  itself  for  area  computation; 
unfortunately,  the  triangulation  may  produce  edges  whose  lengths 
exceed  the  nodes’  communication  radius  by  up  to  a  factor  of  2. 
This  means  that,  while  a  valid  triangle  could  be  discovered  at  one 
vertex,  the  other  pair  of  vertices  in  the  triangle  may  be  out  of  ra¬ 
dio  range  of  each  other  and  hence  will  not  themselves  discover  the 
triangle.  Allowing  a  two-hop  communication  path  between  these 
vertices  does  not  present  a  problem  in  our  application;  indeed, 
paths  between  vertices  will  generally  become  multi-hop  as  the 
transform  scale  becomes  coarser.  Therefore,  we  propose  adding  an 
additional  communication  step  to  the  method  of  [3]:  Once  a  node 
has  identified  its  set  of  valid  triangles,  it  informs  the  other  pair  of 
nodes  in  each  triangle  of  their  membership  in  that  triangle.  Us¬ 
ing  found  and  received  triangles,  a  node  can  then  incorporate  any 
non-included  Gabriel  edges,  forming  additional  triangles.  While 
this  requires  some  additional  negotiation  with  neighbors,  the  num¬ 
ber  of  Gabriel  edges  remaining  after  the  second  round  of  message 
passing  is  typically  very  small  in  practice. 

Ideally,  the  node  density  and  communication  radius  will  both 
be  sufficiently  high  so  that  most  edges  of  the  perfect  DT  have 
lengths  within  twice  that  radius.  When  this  is  the  case,  the  dis¬ 
tributed  triangulation  will  produce  a  mesh  such  as  Figure  1(a), 
where  all  but  the  exceptionally  long  edges  at  the  boundary  of  the 
DT  are  captured  in  the  distributed  mesh.  When  interior  links  of 
the  DT  exceed  this  bound,  the  distributed  mesh  will  exhibit  non- 
triangular  polygons  in  its  interior,  whose  perimeters  create  interior 
“boundaries”  that  must  be  treated  similarly  to  exterior  boundaries 
to  ensure  stability  of  the  predict  operator  of  Section  5.2. 

When  the  node  density  and  communication  radius  allow  the 
distributed  triangulation  to  capture  all  but  the  longest  edges  on  the 
boundary  of  the  DT,  re-triangulation  of  the  mesh  after  decimation 
follows  easily  using  the  technique  of  [15].  Since  only  interior  ver¬ 
tices  are  subject  to  decimation,  such  nodes  can  locally  compute 
new  Delaunay  triangles  and  inform  their  previous  neighbors  of 
their  new  connections  with  each  other,  maintaining  a  Delaunay 
mesh  throughout  scale.  In  the  case  that  the  distributed  triangula¬ 
tion  differs  on  its  interior  front  the  DT,  mesh  refinement  can  pro¬ 
ceed  with  the  predicted  point  linking  its  closest  neighbor  to  the 
others. 


with  the  edge  does  not  contain  any  other  nodes  of  the  triangulation. 


Fig.  5.  Results  of  nonlinear  approximation  experiments.  Mean  square  error  (log-scale)  is  plotted  on  the  y-axis  with  the  number  of  ap¬ 
proximation  coefficients  used  in  the  reconstruction  on  the  x-axis.  We  note  the  rapid  decay  of  the  approximation  error  of  (a)  the  noisy,  dis¬ 
continuous  quadratic  field  and  (b)  the  smooth  Gaussian  bump  field.  Performance  is  comparable  to  more  traditional  biorthogonal  CDF(2,2) 
wavelets  on  square-grid  sampling  for  (c)  the  noisy,  discontinuous  quadratic  and  (d)  smooth  Gaussian  bump  fields. 


8.  DISTRIBUTED  APPROXIMATION  EXAMPLES 

We  present  two  brief  evaluations  showcasing  the  effectiveness  of 
our  proposed  transform. 

First,  echoing  an  example  of  [1],  we  simulate  in-network  non¬ 
linear  approximation  (compression  without  coefficient  quantiza¬ 
tion)  of  the  transform  coefficients  prior  to  transmission  to  an  ex¬ 
ternal  data  sink.  In  this  scenario,  the  sink  floods  the  network  with  a 
magnitude  threshold  below  which  wavelet  coefficients  are  consid¬ 
ered  insignificant.  Each  node  with  a  significant  coefficient  passes 
its  coefficient  value  and  node  ID  to  the  sink,  which  reconstructs  the 
measurement  field  via  the  inverse  wavelet  transform.  We  include 
successively  smaller  coefficients  in  the  reconstruction  to  generate 
nonlinear  approximation  curves  of  increasing  fidelity.  We  con¬ 
sider  the  two  fields  shown  in  Figure  4;  the  first  is  a  noisy  quadratic 
bump  with  a  discontinuity,  and  the  second  is  a  set  of  random 
Gaussian  bumps  populating  a  smoothly  varying  quadratic  field. 
Both  examples  contain  super-planar  features  beyond  the  scope  of 
piecewise-planar  approximation.  For  each  field,  we  created  300 
sample  points  from  uniformly  distributed  random  node  locations. 
Figures  5(a)  and  (b)  display  the  nonlinear  approximation  results 
for  each,  plotting  along  the  y-axis  the  mean-square  approximation 
error  on  a  log  scale  and  along  the  x-axis  the  number  of  wavelet 
coefficients  used  in  the  approximation.  Both  curves  exhibit  the 
smooth,  rapidly  decreasing  decay  we  expect  from  a  wavelet  trans¬ 
form  on  piecewise-smooth  data. 

Second,  as  a  sanity  check,  we  also  compare  the  performance 
of  our  transform  to  a  biorthogonal  wavelet  transform  with  the 
same  number  of  vanishing  moments  (CDF(2,2)  [4])  on  a  regularly 
spaced  square  grid  of  256  grid  points.  We  compare  the  nonlinear 
approximation  curves,  as  shown  in  5(c)  and  (d).  The  solid  lift¬ 
ing  approximation  curve  is  nearly  coincident  with  the  CDF  dashed 
curve,  indicating  that  our  distributed  lifting  transform  is  competi¬ 
tive  with  traditional  wavelet  techniques. 

9.  CONCLUSIONS  AND  FUTURE  WORK 

We  have  developed  a  fully  distributed,  irregular-grid  wavelet  trans¬ 
form  and  protocol  for  sensor  networks  that  is  capable  of  piecewise- 
planar  multiscale  approximation.  We  have  presented  distributed 
solutions  to  implementation  issues  included  mesh  building,  fil¬ 
ter  coefficient  calculation,  and  transform  coefficient  calculation. 
The  transform  coefficients  have  desirable  performance  under 
threshold-based  processing  such  as  distributed  compression  and 
even  perform  competitively  with  centralized,  regular-grid  wavelet 


techniques.  In  current  work,  we  are  developing  a  suite  of  dis¬ 
tributed  processing  algorithms  based  on  this  transform  and  explor¬ 
ing  issues  of  higher  order  approximation,  transform  stability,  trans¬ 
form  coefficient  quantization,  and  extension  to  3-D  grids. 
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